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ABSTRACT 


Increasing  the  nitramine  content  of  solid  rocket  propellants  increases 
the  overall  performance  of  the  system  as  well  as  the  sensitivity  to  detonation 
by  shock  initiation.  In  some  instances  a  confined  zone  of  granulated  propel¬ 
lant  adjacent  to  a  zone  of  cast  propellant  can  provide  a  rapid  enough  pres¬ 
sure-rise  rate  to  shock  initiate  the  cast  material.  If  the  cast  propellant 
is  porous,  the  detonation  will  initiate  at  some  location  ahead  of  the  granu¬ 
lated  bed/cast  material  interface.  The  work  presented  here  is  an  effort 
to  numerically  model  this  Deflagration  to  Shock  to  Detonation  Transition 
(DSDT)  event.  Results  are  presented  showing  the  detonation  build  up  for 
propellants/explosives  with  various  initial  void  content  and  ramp  wave  com¬ 
pression  loads. 


CHAPTER  1 


REVIEW  AND  OVERVIEW 

Improving  the  specific  impulse  of  a  solid  rocket  motor  continues  to  be  an 
engineering  endeavor.  The  specific  impulse  of  a  system  is  defined  to  be  the 
ratio  of  thrust  to  fuel  mass  flow  rate  and  is  a  technical  description  of  the 
overall  performance  of  a  system.  Recently,  one  of  the  most  advantageous  ways 
of  increasing  the  specific  impulse  has  been  to  use  secondary  high  explosives 
as  constituents  in  the  propellant  mixture.  However,  when  high-energy 
nitramines  such  as  HMX  (octogen)  are  utilized  in  the  propellant  formulation 
the  hazard  of  a  Deflagration  to  Shock  to  Detonation  (DSDT)  becomes  a  relevant 
new  concern. 

A  DSDT  event  is  defined  as  a  controlled  subsonic  deflagration  wave  making 
a  transition  to  a  high  order  steady  detonation  wave.  Most  researchers  agree 
that  in  order  for  DSDT  to  occur  the  rocket  motor  grain  must  first  be 
damaged.  Granulation  of  the  propellant  bed  can  be  a  result  of  a  handling 
accident  or  case  and/or  nozzle  failure  during  operation.  Moreover,  once  a 
region  of  granulated  material  has  been  ignited  in  a  confined  configuration  the 
high  surface  to  volume  ratio  particles  provide  an  increased  gas  generation  and 
a  rapid  pressurization  rate  which  can  shock  initiate  the  detonative 
reaction.  Therefore  the  occurrence  of  a  DSDT  event  can  result  in  total 
destruction  of  the  solid  propellant  rocket  motor  assembly  in  the  order  of  only 
milliseconds  after  the  onset  of  granulation. 

Since  the  employment  of  high  secondary  explosives  in  the  propellant 
formulation  is  now  taking  place  (to  increase  the  specific  impulse),  there  is  a 
need  to  understand  the  criteria  for  a  DSDT  event  to  occur  so  the  hazard  can  be 


eliminated  or  at  least  avoided.  In  recent  years  there  has  been  an  increasing 
amount  of  research  in  this  particular  area  of  hazards.  At  the  University  of 
Illinois,  under  the  direction  of  Professor  Herman  Krier,  an  effort  has  been 
made  for  over  a  decade  to  investigate  the  occurrence  of  a  Deflagration  to 
Shock  to  Detonation  Transition  in  cyclotetramethylene  tetranitramine  (HMX, 
octogen).  Examples  of  previous  studies  are  given  in  references  [1-4].  The 
research  in  this  study  is  an  attempt  to  delineate  a  DSDT  event  from  the  onset 
of  compression  waves  propagating  into  the  propellant  bed  to  an  eventual  steady 
detonation. 

1.1  Related  Published  Work 

As  mentioned,  here  at  the  University  of  Illinois,  Krier  and  co-workers 
have  extensively  examined  the  possibility  of  DSDT  occurring  as  a  result  of 
convective  flame  propagation  through  granular  propellants  [1-4],  The 
propagation  sequence  of  events  can  be  summarized  as  follows.  First,  pressure 
gradients  develop  in  the  granular  bed  from  the  localized  burning  of  propellant 
fragments.  To  conserve  momentum,  the  hot  product  gases  are  driven  through  the 
cracks  between  the  propellant  fragments.  Convective  heat  transfer  from  the 
product  gases  to  the  surface  of  particles  ignites  more  fragments,  which  in 
turn  increase  the  pressure  gradients.  As  more  and  more  fragments  are  ignited 
the  pressure  gradients  increase  in  magnitude  eventually  leading  to  shock 
initiation  of  the  detonation  of  the  unreacted  propellant.  This  process  of 
events  is  referred  to  by  some  as  the  accelerated  convective  burn  model, 
developed  by  Butler,  Lembeck  and  Krier  [4],  In  brief,  the  rapid  pressure 
rise,  fueled  by  particle  ignition  from  convective  heat  transfer,  appears  to  be 
the  primary  action  that  leads  to  detonation. 


On  the  other  hand,  Campbell's  investigation  [5]  suggested  that  convective 
flame  spreading  was  essential  only  in  the  early  stages  of  OSDT.  He  suggested 
that  the  stress  waves  produced  by  the  burning  fragments  propagate  ahead  of  the 
flame  front  and  subsequently  cause  detonation.  As  the  pressure  gradients 
increase,  stronger  compression  waves  propagate  and  coalesce  with  earlier  waves 
to  form  and  continually  strengthen  the  shock  prior  to  detonation.  Campbell 
went  on  to  hypothesize  that  once  a  critical  pressure  was  reached  the  pores 
would  collapse  between  the  fragments  and  cause  a  solid  plug  to  form.  The  plug 
would  continue  to  grow  until  the  shock  wave  was  strong  enough  to  initiate 
detonation  and  in  some  instances  retonation. 

To  validate  his  theory  experimentally,  Campbell  packed  a  thick  walled 
steel  pipe  with  granulated  HMX  and  inserted  one  or  more  neoprene  diaphragms  at 
various  locations  throughout  the  interior  of  the  column.  The  neoprene  disks 
completely  covered  the  cross-sectional  area  of  the  pipe,  preventing 
penetration  of  the  hot  product  gases  beyond  the  location  of  the  disk.  For 
various  degrees  of  fineness  of  granulated  HMX,  the  experiments  showed 
detonation  does  occur  ahead  of  the  disk  and  in  advance  of  the  convective  flame 
front.  Figure  1.1  delineates  the  procession  of  the  stress  wave  as  time 
increases,  brought  forth  by  the  burning  propellant  fragments,  in  times  t^ 
through  tg.  At  time  tg,  the  stress  wave  has  advanced  ahead  of  the  first 
disk.  Subsequently,  the  figure  portrays  at  times  tg  and  t^,  the  detonation 
wave  traveling  through  the  granulated  bed.  Furthermore,  for  the  same 
effective  diameters,  experimental  runs  were  made,  absent  of  the  neoprene 
disks,  with  coinciding  run-up  to  detonation  distances.  In  summary,  Figure  1.2 
depicts  the  measured  run-up  length  to  detonation  for  HMX  granulated  beds,  with 
and  without  interval  barriers. 


■t  A:  Schematic  of  Campbell's  experiment  [5]  and  illustration 
assumed  pressure  -  x-location  profiles. 

t  B:  The  burning  granular  propellant  in  Zone  1  shock  initiates 
:  propellant  in  Zone  2. 
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Unlike  Campbell  [5],  Macek  [6]  and  Tarver  [7]  et  al .  studied  DSDT  in 
homogeneous  solid  explosives.  Macek  ran  a  series  of  experiments  in  which 
heavily  confined  cast  cylinders  of  dicthynitramine  dinitrate  (DINA)  and  50/50 
pentolite  were  thermally  ignited  by  an  electrical  source.  The  pressure  rises 
observed  in  the  experiments  by  Macek  were  approximated  by  an  exponential, 
P=0.08  exp(0.1  t),  where  pressure,  P,  and  time,  t,  had  units  of  gigapascals 
and  microseconds,  respectively.  Moreover,  Macek  used  a  simplified  model  to 
serve  as  a  prototype  of  explosive  burning  under  confinement.  The  model 
employed  a  linear  burn  rate  at  the  plane  of  deflagration,  which  is  normal  to 
the  direction  of  the  propagation  of  the  flame,  that  separates  the  product 
gases  from  the  unreacted  explosive.  Accordingly,  the  method  of 
characteristics  was  utilized  to  obtain  theoretical  verification  for  the 
experimental ly  observed  run-up  distance  to  detonation.  The  point  of 
coalescence  of  the  right-running  characteristics,  as  shown  in  Figure  1.3  was 
conjectured  to  be  the  point  of  formation  of  the  shock  which  consequently 
initiates  detonation. 

However,  Jacobs  [8]  found  that  Macek  [6]  had  neglected  to  include  the 
compressibility  of  the  material,  and  with  this  correction  in  the  model,  his 
recent  calculations  infer  that  conductive  burning  could  not  have  produced  a 
shock  wave.  Both  Jacobs  [8]  and  Tarver  C/1  concluded  that  a  mechanical  means 
of  increasing  the  burning  surface  must  exist  to  obtain  exponential  burning 
rates  that  exceed  the  expected  normal  rates  by  several  orders  of  magnitude. 
Hence,  Anderson  and  Kooker  [9]  postulated  that  deconsolidation  of  a  slightly 
porous  material  could  x>ccur  through  confined  burning  by  shear-induced 
stresses,  thus  creating  a  greater  surface  to  volume  ratio. 


Vicinity  of  Shock  Formation 


Figure  1.3  Lines  of  constant  stress  in  homogeneous  material  (HMX)  being 
stressed  at  x=0,  taken  from  Reference  [10]. 


Recently  Coyne,  Butler  and  Krier  [10]  studied  the  propagation  of  stress 
waves  into  porous  HMX.  During  their  investigation  they  found  it  difficult  to 
employ  the  method  of  characteristics  while  utilizing  a  version  of  the  more 
realistic  Mie-Grunsien  equation  of  state.  Instead  they  used  the  Tait  equation 
of  state,  which  does  not  properly  represent  the  isentrope,  to  obtain  a 
numerical  solution  of  the  conservation  equations  by  a  Lagrangian  finite 
differencing  technique.  This  solution  was  to  be  verified  by  comparing  results 
to  those  obtained  by  the  method  of  characteristics.  Depicted  in  Figure  1.4 
are  the  results  of  the  comparison.  Eventually  the  modified  Mie-Gruneisen 
relation  was  incorporated  in  the  finite  difference  code  to  model  stress  wave 
propagation  in  a  porous  nonreactive  material. 

The  "Pop-plot"  is  named  after  its  originator,  Nickalous  Popaloto,  and 
delineates  the  shock  pressure  strength  to  run-up  distance  to  detonation  on  a 
log-log  scale.  Figure  1.5  is  an  example  of  a  "Pop-plot"  which  was  obtained  by 
Dick  [11]  for  porous  HMX.  Due  to  the  hazards  associated  when  experimenting 
with  explosives  and  the  difficulties  in  obtaining  accurate  results,  there  is 
little  data  on  run-up  distances  to  detonation  for  granular  explosives. 

However,  using  a  data  acquisition  technique  different  from  the  standardized 
wedge  test,  J.  J.  Dick  obtained  "Pop-plot"  data  for  porous  samples  of  HMX  of 
initial  density  pQ=  1.24  +  D.04  g/cm^  [111.  The  wedge  test  records  the 
trajectory  of  a  shock  as  the  wave  travels  through  a  wedged  shaped  sample  by 
photographic  techniques.  Different  from  the  standardized  test,  Dick  measured 
only  the  total  time  that  the  wave  resided  in  the  cylindrical  sample,  thereby 
producing  a  "Pop-plot"  through  extensive  runs  and  statistical  analysis. 

Setchell  [12]  also  ran  experiments  on  initiation  behavior  of  granular 
explosives.  Utilizing  Laser  Velocity  Interferometry,  Setchell  studied  "ramp" 
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Comparison  of  Stress  Wave  Propagation 
using  the  Tait  E.O.S. 


- Finite  Difference  Calculation 

- Method  of  Characteristics 

P[x0(t),t]  =  [0.008  GPa]  exp[(0.1/£is)t] 


90  /msec 


XLOC  (CM) 

Figure  1.4  Stress  wave  propagation  through  solid  HMX,  from  both  the 
finite  difference  calculation  and  the  method  of 
characteristics.  Both  utilizing  the  modified  Tait 
equation  of  state.  Figure  taken  from  Reference  [10]. 


Input  Stress  (GPa) 


Figure  1.5  Run-up  distance  to  detonation  versus  input  stress 
for  several  high-explosives.  Uncertainty  bars  for 
run-up  distance  represent  one  and  two  standard 
deviations.  Data  for  HMX  (1.24  g/cc)  is  from 
Reference  [11].  Data  for  PETN,  PBX-9404,  HMX  (1.89  g/cc) 
is  from  Reference  [31]. 


"ramp"  waves  having  a  finite  rise-time  of  either  0.3  or  0.8  microseconds.  By 
shock  loading  proyoceram,  a  material  known  to  have  stress-strain  relation  with 
negative  curvature,  he  was  able  to  disperse  a  shock  wave  into  a  wave  having  a 
finite  rise  time.  Setchell  found,  by  comparing  velocity-time  profiles  of 
tests  on  P8X-9404  shocked  to  the  same  peak  pressure,  that  very  little 
chemical  energy  was  released  in  the  tests  with  ramp  waves  prior  to  shock 
formation.  His  records  indicate  that  local  hotspot  generation  and  ignition 
are  strongly  inhibited  by  finite  rates  of  compression. 

1.2  Hot  Spot  Theory 

Although  the  theory  of  hot  spots  as  a  source  for  ignition  in  porous 
material  is  widely  accepted,  there  is  very  little  conclusive  evidence  on  the 
matter  of  generation.  Some  of  the  hypotheses  for  hot  spot  generation  are 
shear  banding,  jetting,  shock  focusing,  and/or  adiabatic  compression.  For 
more  detailed  information,  see  to  References  [13-16].  Hayes  [17]  suggests 
that  the  total  energy  deposited  by  the  shock  wave  can  be  equated  on  a  mass 
fraction  basis  to  the  sum  of  the  reversible  work  done  in  isentropically 
compressing  the  bulk  of  the  material  plus  the  irreversible  heating  of 
localized  hot  spots,  i.e., 

P  +  P 

— 2 -  (vT0-vt)  =  ^H  e^vT*^H^  +  (1*1) 

In  equation  (1.1)  the  left-hand  side  represents  the  total  energy  deposited  in 
the  material  by  the  shock  of  strength,  P.  The  term  e j ( P )  represents  the 
energy  required  to  isentropically  compress  the  bulk  of  the  material  to  the 
final  shock  pressure,  and  the  remaining  energy  term,  e(vT,TH),  is  the  energy 


available  to  irreversibly  heat  the  hot  spots.  The  Hayes  model  assumes  the 
mass  fraction  of  the  hot  spots,  WH,  to  be  equal  to  the  preshock  specific 
volume  fraction  of  pores 


WH  =  vT0^vS0  "  1  (1<2) 

Here,  the  subscript  'TO'  represents  the  initial  porous  state,  and  the 
subscript  'SO'  refers  to  the  homogenous  initial  state.  Furthermore,  Hayes 
[17]  states  that  scissing  of  particular  chemical  bonds,  stemming  from  the  high 
frequency  content  in  the  shock  front  and  the  construction  of  a  non-equilbrium 
temperature,  can  be  responsible  for  the  increase  in  reactivity  and 
subsequently  the  decrease  in  decomposition  times  for  the  hot  spots  observed  in 
his  experiments.  Figure  1.6  compares  the  observed  decomposition  times  for 
hexanitrostilbene  (HNS)  to  times  derived  from  low-temperature  Arrhenius 
kinetics  for  a  hot  spot  temperature  regime.  Consequently,  the  explosion  model 
fails  to  predict  corresponding  decomposition  times  and  falls  away  by  a  bigger 
margin  as  temperature  of  the  hot  spots  decrease. 

1.3  Scope  of  Our  Study 

As  previously  mentioned,  only  if  there  is  fragmentation  of  the  propellant 
bed  in  the  burning  region,  along  with  proper  confinement,  will  there  be 
sufficient  pressurization  rates  to  shock  initiate  detonation  for  cast 
explosives  [18].  Consider  a  rocket  motor  with  a  center  burning  grain 
configuration  as  shown  in  Figure  1.7.  To  illustrate  a  DSDT  event,  a  crack  in 
the  propellant  grain  is  assumed,  and  within  the  crack  a  packed  bed  of 
fragments  of  known  surface-to-volume  ratio  exist.  The  fragment  filled  crack 
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Fiqure  1.6  Experimentally  observed  and  calculated  hot  spot 

decomposition  times  versus  the  hot  spot  temperature 
inverse  for  HNS.  Figure  taken  from  Reference  [l/J. 


DSD T  HAZARD  IN  ROCKET  MOTOR 

Figure  1.7  Sketch  of  solid  propellant  rocket 
motor  with  a  crack  indicated  in  the 


is  enlarged  and  shown  in  Figure  1.8.  At  some  arbitrary  time,  the  bed  ignites 
and  starts  burning  at  x=0  (left-end  of  figure).  Increased  product  gas 
generation  beyond  the  level  necessary  for  steady  state  motor  operation 
orginates  from  the  assumed  high  surface-to-vol ume  ratio  of  confined  sub¬ 
millimeter  size  particles.  In  other  words,  the  amount  of  gas  being  generated 
by  the  decomposing  propellant  far  exceeds  the  amount  exiting  the  nozzle. 

If  the  length  of  the  packed  bed  is  longer  than  the  critical  condition  for 
accelerated  convective  combustion  to  occur,  subsequent  detonation  is 
inevitable  [19].  However,  if  the  bed  length  is  less  than  the  critical  length, 
pressure  gradients  produced  by  the  granulated  propellant  can  provide  the 
impetus  to  shock  initiate  detonation  in  the  adjacent  region  of  cast  explosive, 
depicted  as  zone  1  in  Figure  1.9.  Only  stress  waves  can  be  transmitted  across 
the  zone2/zonel  interface.  Even  though  the  solid  may  contain  pores,  it  is 
assumed  to  be  impermeable  to  the  flow  of  hot  gases.  A  schematic 
representation  of  the  sequence  of  events  leading  to  Deflagration  to  Shock  to 
Detonation  Transition  modeled  in  this  study  is  shown  in  Figure  1.10.  A  value 
of  <p  equal  to  unity  represents  a  zone  of  all  gas,  while  $  equal  to  zero 
indicates  a  homogeneous  solid.  In  Figure  1.10,  the  heavy  black  dots  are 
representative  of  microvoids  in  the  cast  material.  Illustrated  in  Part  B  is 
the  collapse  of  the  pores,  a  result  of  the  stress  load  transmitted  across  the 
granulated  bed/cast  explosive  interface.  Parts  C  and  D  show  the  length  of  the 
pore  collapse  zone  to  increase  with  time  as  the  lead  compression  waves  travel 
farther  into  the  explosive.  The  finite  compression  waves  coalesce  into  a 
shock  front  which  then  initiates  the  cast  explosive  downstream  of  the 
interface.  From  this  location  a  detonation  wave  propagates  through  the 
porous  material,  while  a  retonation  wave  propagates  back  through  the 
compressed  material  (Part  E). 


The  purpose  of  the  research  study  here  is  then  to  model  the  key  elements 
of  the  five  part  scenario.  A  one-dimensional  hydrodynamic  Lagrangian  finite 
difference  technique  is  used  to  numerically  solve  the  conservation  equations 
of  mass,  momentum,  and  energy.  A  pore  collapse  theory  derived  by  Carroll  and 
Holt  [20]  which  demarcates  three  regimes  of  deformation,  elastic,  elastic- 
plastic,  and  plastic  is  utilized  to  determine  the  rate  of  compaction  and  the 
development  of  the  solid  plug.  Rasically,  this  portion  of  the  code  is  an 
extension  on  previous  work  done  by  Coyne  [21].  Furthermore,  the  Hayes  hot 
spot  theory  is  incorporated  in  the  code  to  define  the  sensitivity  to 
reaction.  By  introducing  reactive  chemistry  to  the  code,  a  strong  effort  is 
made  to  model  the  detonation  and  retonation  waves  which  are  initiated  by  a 
shock  wave  generated  from  ramp  wave  inputs  with  rise  times  on  the  order  of 
tens  of  microseconds.  It  was  shown  in  Reference  [19]  that  rise  times  of  this 
order  are  typical  for  burning,  granulated  beds  whic1’  have  lengths  less  than 
their  critical  detonation  run-up  length.  Although,  Setchell  [12]  also  studied 
ramp  waves,  the  rise  times  were  faster  by  an  order  of  magnitude. 
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CHAPTER  2 


MATHEMATICAL  AND  NUMERICAL  ANALYSIS 


The  equations  which  mathematically  model  shock  initiation  to  detonation 
due  to  a  rapid  compression  of  high  explosive  containing  voids  are  presented  in 
the  sections  to  follow.  Specifically,  the  material  being  considered  here  is 
HMX,  requiring  specific  constitutive  relations.  Properties  for  HMX  are  listed 
in  Appendix  A.  The  model  clearly  can  be  applied  to  other  reactive  solids,  but 
an  equation  of  state,  a  caloric  relation,  and  the  material  Hugoniot  must  be 
known. 


2.1  Free  Boundary 

The  first  step  taken  to  model  DSDT  in  a  porous  explosive  solid  was  to 
determine  the  pressure  gradients  produced  by  the  adjacent  reacting  granulated 
bed.  Zone  2  in  Figure  1.9.  An  analysis  by  Butler  and  Krier  described  in 
Reference  [4],  models  the  accelerating  convective  burn  and  rapid  pressurizaton 
for  granulated  beds  of  explosives. 

Figure  2.1  presents  the  predicted  pressure  rise  rate  at  one  location  in 
such  a  bed.  The  bed  length  is  always  less  than  the  required  run-up  distance 
to  detonation,  £Cj.  Notice  that  the  pressurization  rate  is  strongly  dependent 
on  particle  size,  increasing,  as  one  would  expect,  with  smaller  propellant 
fragments.  Subsequently,  the  rate  at  which  the  impermeable  bed  is  being 
stressed  (assuming  that  such  a  bed  is  adjacent  to  the  porous,  permeable  bed) 
is  defined  by  the  pressure-time  functions  predicted  by  the  Deflagration  to 
Detonation  Transition  case  documented  in  Reference  [19]. 


1 


The  pressure-time  function  here  is  approximated  by  a  linear  p-t  relation. 


with  the  slope  a  specific  function  of  the  material  properties.  Therefore  to 
simplify  the  loading  boundary  condition  at  the  granulated/cast  bed  explosive 
interface,  the  pressure  magnitude  of  the  left  free  boundary  was  assumed  to 
satisfy 


P(t)  -  (P*  -  P0)(t/t*)  +  P0  t<_t*  (2.1a) 

P(t)  =  P*  t  >_  t*  (2.1b) 

Here  P*  is  the  maximum  pressure  applied  to  the  left  boundary,  while  t*  is  the 
characteristic  rise-time  for  the  "ramp"  wave  to  reach  the  maximum  stress  at 
the  free  boundary.  It  should  be  noted  that  a  typical  range  of  t*  includes  the 
interval  1  <  t*  <  50  microseconds.  Thus  this  type  of  "ramp-loading"  is  much 
slower  than  the  "ramp"  wave  compressions  of  explosives  studied  by  Setchell 
[12],  where  0.3  <  t*  <  0.8  microseconds.  A  shock  loading  to  P*  is  usually 
assumed  to  be  t*  <  0.01  microseconds. 

2.2  Governing  Equations 

For  our  hydrodynamic  analysis  the  Lagrangian  or  material  form  of  the 
governing  equations  was  chosen,  instead  of  Eulerian  form.  The  Lagrangian 
coordinates  are  fixed  to  the  material  and  follow  this  material  as  it  moves 
with  time,  whereas  the  Eulerian  coordinates  are  a  fixed  frame  of  reference, 
where  mass,  momentum,  and  energy  may  enter  or  exit  through  the  control 
surfaces.  Therefore,  with  the  Eulerian  formulation,  the  boundary  location 
must  be  implicitly  determined  after  each  time  increment.  However,  with  the 
Lagrangian  description,  the  location  of  the  free  boundary  condition  is 
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explicitly  known,  since  the  Lagrangian  boundary  coordinates  move  with  the 
material  boundary. 

The  inviscid  Lagrangian  one-dimensional  unsteady  form  of  the  conservation 
of  mass,  momentum,  and  energy  equations  are  expressed  for  total  mechanical 
mixture,  respectively  as 


pll+o 

at  w 


(2.2) 

(2.3) 

(2.4) 


In  equation  (2.2  -  2.4),  v  represents  the  specific  volume;  u,  particle 
velocity;  e,  the  specific  internal  energy;  P,  the  total  stress;  and  Q,  the 
heat  added  by  chemical  reaction  per  unit  mass  per  unit  time.  The  symbols 
JR  and  indicate  the  partial  derivatives  with  respect  to  the  Lagrangian 
spatial  coordinate  and  time,  respectively.  For  comparison  purposes,  the 
Euleri an  form  of  the  conservation  of  mass,  momentum,  and  energy  equation, 
respectively,  would  be. 
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where  the  symbol  —  denotes  the  partial  derivative  with  respect  to  the 

a  X 

Eulerian  spatial  coordinate.  Notice  that  equations  (2.5  -  2.7),  in  absence  of 
the  convective  flux  terms,  are  identical  to  equations  (2.2  -  2.4).  Since  the 
same  mass  is  contained  in  a  given  Lagrangian  control  volume  throughout  the 
variation  of  time  the  convective  flux  terms  do  not  appear  in  equations  (2.2  - 
2.4). 

Generally,  the  Lagrangian  form  of  the  governing  equations  are  easier  to 
incorporate  into  a  numerical  solution  technique,  since  it  is  clear  that  they 
are  in  a  simpler  form  than  the  Eulerian,  and,  most  importantly,  one  has  the 
advantage  when  dealing  with  a  free  or  moving  boundary  condition.  To  provide 
mathematical  closure  to  the  governing  equations,  for  an  inert  material,  an 
equation  of  state  is  needed,  since  the  three  equations  (2.2)  to  (2.4)  involve 
four  unknowns. 


2.3  Constitutive  Relations 

From  the  Second  Law  of  Thermodynamics,  there  exists  a  unique  relationship 
between  the  material  equation  of  state,  Ps(v$,  Ts),  and  the  caloric  equation 
of  state,  es(vs,  T$).  One  way  of  relating  the  two  state  equations  is  through 
Helmholtz  free  energy  function  and  its  thermodynamic  derivatives.  Therefore 
Ps(vs,  Ts)  and  es(v$,  T$)  must  satisfy  the  reciprocity  relations. 


P 


s 


(2.8) 


t  ii_ 
Ts  3TS 


(2.9) 


In  equations  (2.8)  and  (2.9),  T  represents  temperature,  P,  pressure,  e, 

specific  internal  energy,  v,  specific  volume,  and  \j>  Helmholtz  free  energy. 

The  subscript  s  denotes  a  solid  material,  while  the  symbols 
3  3 

—  and  jy  indicate  partial  derivatives  with  respect  to  specific  volume  and 
temperature,  respectively.  Appendix  B  gives  a  review  of  these  fundamental 
thermodynamic  concepts. 

With  the  assumption  that  the  Gruneisen  coefficient,  r,  is  constant,  the 
Helmholtz  free  energy  function  is  expressed  in  the  following  form  (Baer  and 
Nunziato  [22]). 


where  cvs  represents  the  specific  heat  at  constant  volume  of  the  solid  phase 
which  is  also  assumed  to  be  a  constant.  The  term  J(v$)  in  equation  (2.10)  is 
a  nonlinear  volume-dependent  function  determined  from  shock  Hugoniot 
experiments  [22]. 

The  Gruneisen  coefficient  is  defined  by  the  thermodynamic  derivative 

r  (v)  =  -  v  (|£)v  (2.11) 

which  characterizes  the  ratio  of  thermal  pressure  to  the  thermal  energy  of  the 
lattice.  At  standard  volume,  the  Gruneisen  coefficient,  r0  s  r(vQ),  of  a 
material  can  be  related  to  other  properties,  such  as  the  isothermal 


compressibility  of  the  material  at  standard  conditions,  tco 
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_i  (11) 

vQ  aP  t 


(2.1?) 


and  the  coefficient  of  thermal  expansion  at  constant  pressure,  <p 


1  fill 

*P  '  vQ  3T  p 


(2.13) 


Thus  an  expression  can  be  obtained  for  the  Gruneisen  coefficient  at 
standard  volume,  i.e.. 


v  K 
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(2.14) 


where  c0  is  the  ambient  sound  speed  of  the  material. 

Although  equations  (2. 2-2. 4,  2.8  and  2.9)  form  a  basis  for  the 
mathematical  description  of  a  homogeneous  material,  an  additional  relation  is 
needed  to  describe  a  heterogeneous  material.  Porosity  is  defined  as  the  ratio 
of  total  volume  of  the  material  to  the  volume  occupied  by  the  solid  phase,  or 


(2.15) 


Initially,  the  volume  not  occupied  by  the  solid  material  is  assumed  to  be 
massless.  With  the  introduction  of  porosity  an  additional  equation  is  needed 
to  complete  the  mathematical  description. 

Extensive  research  in  the  area  of  mathematically  modeling  the  collapse  of 
a  porous  material  under  an  applied  external  load  has  been  performed  by  Carroll 
and  Holt  [20].  In  their  model  the  porous  matrix  is  exemplified  as  a  hollow 


sphere  where  the  inner  and  outer  radii  are  chosen  such  that  the  overall 
porosity  of  the  material  is  accurately  portrayed.  As  previously  mentioned  in 
Chapter  1,  the  model  of  Carroll  and  Holt  [20]  assumed  pore  collapse  to  occur 
in  three  distinct  regimes:  (1)  elastic  phase,  where  elastic  deformation  in 
the  solid  takes  precedence,  (2)  elastic-plastic  phase,  where  plastic 
deformation  initially  starts  at  the  inner  radius  and  subsequently  progresses 
outward  until  plastic  deformation  begins  to  occur  at  the  outer  radius,  and 
(3)  plastic  phase,  where  plastic  deformation  occurs  throughout  the  sphere. 
The  so  called  "P  -  a  relations"  for  the  three  particular  phases  of  compaction 
and  the  appropriate  range  over  which  each  applies  are  given  by. 


elastic  phase  «o  2.  a  2.  ai 

4  G(a0  -  a) 

P  3  a (a  -  1) 


elastic-plastic  phase  aj  ®  2.  ap 


P  =  7  Y  t1  '  H  (<V  a)  +  *n["Y(a°- 


2G(a  -  a) 


plastic  phase  ag  >_  a  >  1 


P  =  |  Y  *n  (^-fy) 


(2.16a) 


(2.16b) 


(2.16c) 


where  the  limits  between  the  three  phases  are  expressed  as 
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a2  =  2G  +  Y 


(2.18) 


Kooker  and  Anderson  [9]  have  also  used  the  static  pore  collapse  model  with  the 
yield  stress,  Y,  and  shear  modulus,  G,  expressed  as  functions  of  a,  namely. 


Y  =  Yq/(2c x  -  1)‘ 


(2.19) 


G  =  G0  exp  [-5(o  -  l)/a] 


( 2.20) 


Here,  Y0  and  G0  are  the  initial  yield  stress  and  shear  modulus, 
respectively.  Figure  2.2,  taken  from  Reference  [21],  illustrates  the  three 
regimes  of  deformation  for  various  initial  porosities.  Notice  that  the 
plastic  phase  approaches  unity,  complete  compaction.  During  deforma tive 
compression  the  pressure  of  the  mechanical  mixture  s  equated  to  the  pressure 
of  the  solid 


P  =  Pc 


(2.21) 


Thus,  with  the  additional  parameter  a  and  the  "P  -  a  law",  a  mathematical 
description  for  a  nonreactive  porous  material  is  complete. 

The  governing  equations  for  the  conservation  of  mass,  momentum,  and 
energy  represented  in  Section  2.2  are  written  in  terms  of  the  thermodynamic 
properties  (P,v,T,e)  of  the  mechanical  mixture  as  well  as  the  dynamic  variable 
u.  However,  during  reaction  additional  relations  are  needed  to  separate  the 
individual  phase  properties  of  the  solid  and  product  gases  from  those  of  the 
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Figure  2.2  Pressure-void  volume  relationship  for  porous  HMX. 

Figure  taken  from  Reference  [21]. 
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mixture.  Hence,  it  is  conjectured  that  any  arbitrary  volume,  V y,  within  the 
continuum  can  contain  both  the  solid  and  gas  phases.  Therefore  the  individual 
phase  volumes  sum  to  the  total  volume,  i.e. 


VT  =  Vg  +  V$  (2.22) 

Analogously,  the  total  energy  of  the  mechanical  mixture,  Ey,  is  the  sum  of  the 
total  energy  of  both  the  solid  and  gas  phase,  represented  as 


i 


ET  -  Eg  *  Es  (2.?3) 

With  the  introduction  of  a  product  gas  phase  in  the  mechanical  system, 
another  constitutive  equation  must  be  prescribed.  A  nonideal  covolume 
equation  of  state  was  chosen,  similar  to  the  one  utilized  by  Butler  and  Krier 
[19],  for  the  product  gas  phase,  namely 


Pg  ■  iTg  (1  f  S/»g)/V9 


(2.24) 


where  R  is  the  gas  constant,  and  g  is  a  covolume  correction  term.  The  value 
of  g  is  determined  from  values  for  pressure,  temperature,  and  density  at  the 
Chapman-.Jouguet  (C.1)  state  predicted  by  a  thermochemical  code,  TIGER  [23], 
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A  listing  of  the  C.)  values  for  several  loading  densities  of  HMX  is  given  in 
Table  [2.1].  In  accordance  with  the  reciprocity  relations  defined  earlier 
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TABLE  2.1 

CJ  PARAMETERS  PREDICTED  BY  TIGER  [23] 


a0 

PTO 

PCJ 

TCJ 

VCJ 

D 

6 

(g/cc) 

(GPa) 

(K) 

(cc/g) 

(mrn/ys) 

(cc/g) 

1 

1.9 

35.99 

3714 

0.4073 

9.154 

5.192 

1.056 

1.8 

31.70 

3833 

0.4266 

8.711 

4.815 

1.118 

1.7 

27.97 

3931 

0.4478 

8.301 

4.522 

1.188 

1.6 

24.72 

4009 

0.4713 

7.925 

4.300 

1.267 

1.5 

21.86 

4067 

0.4977 

7.581 

4.139 

1.357 

1.4 

19.30 

4106 

0.5278 

7.267 

4.034 

1.462 

1.3 

17. QO 

4126 

0.5627 

6.979 

3.982 

(Equations  2. 8-2. 9),  the  caloric  equation  of  state  for  the  product  gas  phase 
is  expressed  as 

eg  =  cvg  "  ^go^  (2.26) 

where  cVg  is  an  assumed  constant  which  represents  the  specific  heat  at 
constant  volume  of  the  gas  phase.  Furthermore,  while  chemical  reaction  is 
present  the  additional  assumptions, 

P  =  Ps  »  Pg  (2.27) 

T  =  Ts  =  Tg  (2.28) 

are  imposed. 


2.4  The  Local i zee  Hot  Spot  Temperature 

Before  specifically  discussing  the  hot  spot  temperature,  it  is  first 
important  to  compare  the  relative  amount  of  energy  associated  with  shocking  a 
porous  material  compared  to  that  associated  with  shocking  a  homogeneous 
material.  The  text  by  Zeldovich  and  Raizer  [24]  provides  detailed 
background.  Figure  2.3  taken  from  Reference  [21]  delineates  the  Hugoniot  for 
both  a  porous  and  homogeneous  matrix  of  HMX  compressed  to  a  volume  ratio  of 
v/v$0  =  0.9.  In  Figure  2.3  the  horizontally  shaded  area  ABC  and  the 
crosshatched  area  A'B'C  show  the  increase  in  energy  of  a  shocked  porous  and 
homogeneous  material,  respectively.  The  significant  increase  in  energy 
associated  with  shocking  a  porous  material  compared  to  a  homogeneous  material 
is  conjectured  to  be  the  cause  for  a  porous  material  to  be  more  sensitive  to 
shock  initiation  of  the  detonative  reaction  than  a  homogeneous  material. 
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Referring  back  to  Figure  1.5,  a  "Pop-plot"  for  various  propellants,  one  can 
see  verification  of  this.  Note  the  significant  decrease  in  run-up  distance  to 
detonation  for  a  porous  material  compared  to  a  nearly  homogeneous  material, 
shocked  to  the  same  peak  pressure.  In  additon,  triangle  ABA',  shown  in  Figure 
2.3,  is  representative  of  the  amount  of  irreversible  energy  deposited  by 
shocking  a  porous  material. 

As  mentioned  previously  in  Chapter  1,  the  theory  developed  by  Hayes  [17] 
suggests  that  the  total  energy  deposited  by  the  shock  can  be  equated  by  the 
reversible  work  done  in  compressing  the  bulk  material  plus  the  irreversible 
heating  of  the  localized  hot  spots.  The  Equations  (1.1  and  1.2)  underlying 
the  theory,  are  rewritten  as 

P+P 

~2  ^ vTo”  VT ^  =  WH  es^Vs*TH^  +  ^“WH^es^Vs,TI ^  (2.29) 

y 

and  WH  »  -12-  -  1  (2.30) 

H  vso 

where  TH  represents  the  hot  spot  temperature.  For  a  porous  material,  the 
total  energy  of  the  mechanical  mixture  is  assumed  to  equal  the  total  energy  of 
the  solid.  Therefore,  the  isentrope  energy,  ej(P),  can  be  expressed  as 

ej(P)  =  es(v$,Tj)  (2.31) 

where  Tj  is  the  temperature  defined  by  isentropic  compression.  A  relation  for 
the  isentropic  temperature  is  obtained  from  Reference  [24]  and  is  written  as 


With  this  formulation  one  can  calculate  the  hot  spot  temperature.  A  model  is 
now  needed  to  incorporate  this  hot  spot  temperature  into  the  kinetics  which 
represent  how  rapidly  the  solid  explosive  will  gasify  (explode). 

2.5  The  Combustion  Model 

Following  the  simplification  used  by  Mader  [25],  a  first  order  Arrhenius 
burn  model  was  chosen  to  describe  the  chemical  reaction  (decomposition)  rate 
of  the  solid  material,  expressed  as 

~  =  -  z  W  exp  (-EVRT*)  (2.33) 

In  Equation  (2.33)  W  denotes  the  mass  fraction  of  unreacted  explosive,  z, 
frequency  factor,  E*.  activation  energy,  W  ,  universal  gas  constant,  and  T* 
the  characteristic  burn  temperature.  The  symbol  ^  indicates  a  total 
derivative  with  respect  to  time. 

At  first  the  characteristic  burn  temperature  was  represented  by  the 
localized  hot  spot  temperature,  but  Equation  (2.33),  utilizing  values  of  z  and 
E*/  R-  (referred  to  as  activation  temperature)  obtained  from  Nunziato  [26], 
resulted  in  too  slow  of  a  reaction  rate.  Hayes  [17]  also  came  to  this 
conclusion,  based  on  experiments  with  a  similar  propellant,  HNS.  Since  no 
decomposition  times  for  hot  spots  are  available  for  HMX,  we  therefore  fit  our 
model  directly  to  Hayes  data.  By  extracting  several  points  from  the  curve  in 
Figure  1.6,  decomposition  time  was  found  to  be  a  parabolic  function  of  hot 
spot  temperature,  i.e.. 


log  t  =  -0.6744703  2.2482343  [^P-]  -  1.9132332  (2.34) 

'H  H 

where  t  represents  the  decomposition  time  in  microseconds.  The  data  used  by 
Hayes  showed  that  a  detonation  occurred  for  a  corresponding  temperature  of  619 
K  but  did  not  for  a  corresponding  lower  temperature  of  561  K,  therefore  600  K 
was  the  temperature  selected  for  the  lower  limit  in  applying  expression  (2.34) 
in  the  decomposition  model. 

To  incorporate  the  decomposition  time  determined  by  Equation  (2.34)  into 
the  combustion  model.  Equation  (2.33)  was  integrated  to  define  a  new 
activation  temperature  corresponding  to  the  predicted  hot  spot  decomposition 
time.  Therefore,  by  integrating  Equation  (2.33),  the  following  expression  is 
obtained 


(2.35) 


where  Wg  is  the  mass  fraction  remaining  after  chemical  decomposition  of  the 
hot  spot  has  occurred.  8earing  in  mind  that  this  model  is  based  on  HNS  rather 
than  HMX,  the  activation  temperature  predicted  from  Equation  (2.35)  is  bounded 
by  the  value  given  by  Nunziato  [26].  Subsequently,  after  the  hot  spot  has 
decomposed,  the  characteristic  burn  temperature  is  defined  by  one  of  the 
following  two  values;  (1)  the  bulk  temperature  or  (2)  the  average  of  the 
bulk  and  hot  spot  temperature.  These  are  then  two  distinct  cases  which  are 
referred  to  as  model  CB1  and  CB2,  and  are  expressed  below  (along  with  the 
limits  in  which  they  apply)  as 
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Model  CB1 


WB  1  W  >_  ° 


T*  =  T 


(2.36) 


Model  C82 


wB  >_  W  >_  0 


T*  =  (Th  +  T)/2.0  if  TH  >  T 
T*  =  T  if  TH£T 


(2.37a) 

(2.37b) 


Both  models  were  used  in  the  numerical  solution  technique.  Based  on  the 
interpretation  of  the  predicted  results,  model  CB2  was  eventually  assumed 
better.  A  comparison  of  such  results  are  presented  in  the  next  chapter  along 
with  the  explanation  of  why  model  CB2  was  chosen  over  C81. 


2.6  Review  of  the  Key  Assumptions 

Before  proceeding  onto  the  description  of  the  numerical  solution 
technique  a  list  of  the  key  assumptions  made  in  each  portion  of  the  solution 
technique  is  presented  below. 


Homogeneous  Material 


1) 

2) 


Solid  Phase 

3P 

the  Gruneisen  coefficient,  r  =  -  v  (— )y  , 
the  specific  heat  at  constant  volume,  cy$. 


is 


constant. 


is  constant 


Gas  Phase 


1)  the  specific  heat  at  constant  volume,  cVg,  is  constant. 

Porous  Material 

3P 

1)  the  Gruneisen  coefficient,  r  =  -  v  (— )  ,  is  constant. 

2)  the  specific  heat  at  constant  volume,  cys,  is  constant. 

3)  the  voids  initially  in  the  material  are  regarded  as  being  massless. 

Reactive  Material 

1)  the  Gruneisen  coefficient,  r  =  -  v  (-—■)  ,  is  constant. 

d  G  V 

2)  the  specific  heat  at  constant  volume,  cys,  is  constant. 

3)  the  specific  heat  at  constant  volume,  cvg,  is  constant. 

4)  reaction  does  not  initiate  until  the  voids  initially  within  the  matrix  of 
the  material  have  completely  collapsed. 

5)  the  temperature  of  the  gas  phase  equals  the  temperature  of  the  solid  phase 
which  equals  the  temperature  of  the  mechanical  mixture 

6)  the  pressure  of  the  gas  phase  equals  the  pressure  of  the  solid  phase  which 
equals  the  pressure  of  the  mechanical  mixture 

Consequently,  the  subscripts  of  the  parameters  just  equated  will  be  dropped  in 
the  section  which  follow. 

Table  [2.2]  presents  a  method  of  procedure  for  solving  the  flow  equations 
which  summarizes  the  logic,  the  manner  in  which  these  assumptions  constrain 
the  results,  and  the  iteration  pathways  to  assure  convergence  of  all  dependent 
variables. 


TABLE  2.2 


PROCEDURE  FOR  SOLVING  THE  FLOW  EQUATIONS 

1)  Solve  for  the  conservation  of  mass,  Equation  (2.2). 

2)  Solve  for  the  conservation  of  momentum.  Equation  (2.3). 

3)  Solve  for  W,  the  mass  fraction. 

(i)  If  the  pores  initially  contained  in  the  matrix  of  the  material  have 
not  collapsed,  W=1  (no  reaction). 

(ii)  If  the  initial  pores  have  collapsed  and  the  hot  spot  temperature, 
calculated  from  Equation  (2.29),  is  less  than  600  K,  use  Equation 
(2.33)  to  obtain  a  value  for  W. 

(iii)  If  the  initial  pores  have  collapsed  and  the  hot  spot  temperature  is 
greater  than  600  K,  use  Equations  (2.33),  (2.34),  and  (2.35)  to 
determine  a  value  for  W. 

4)  Solve  for  the  conservation  of  energy.  Equation  (2.4). 

(i)  Homogeneous  Material 

(1)  Solid  Phase  -  Solve  for  the  remaining  thermodynamic  parameters, 
(P,T),  knowing  that  the  total  internal  energy  of  the  mechanical 
mixture  is  equal  to  the  total  internal  energy  of  the  solid, 
while  utilizing  the  assumption  that  r  and  cv$  are  constants. 

(2)  Gas  Phase  -  Solve  for  the  remaining  thermodynamic  parameters, 
(P,T),  knowing  the  total  internal  energy  of  the  mechanical 
mixture  equals  the  total  internal  energy  of  the  gas,  while 
making  use  of  the  assumption  that  c„„  is  a  constant. 


(ii)  Porous  Material  -  Solve  for  the  remaining  thermodynamic  parameters, 
(p»vs,T).  Utilizing  the  assumption  that  the  initial  voids  are 
massless  implies  that  the  total  internal  energy  of  the  mechanical 
mixture  is  equal  to  the  total  internal  energy  of  the  solid. 
Furthermore,  an  iterative  technique  must  be  used  to  equate  the 
pressure  predicted  by  the  "P-a"  law  to  the  pressure  of  the  solid  by 
varying  the  specific  volume  of  solid. 

(iii)  Reactive  Material  -  Solve  for  the  remaining  thermodynamic  parameters, 
(P,Vg,vs,T),  by  employing  the  assumption  that  the  temperature  of  the 
gas  and  solid  are  equal  to  each  other.  In  addition,  the  pressure  of 
the  gas  and  solid  can  be  equated  by  utilizing  an  iterative  method 
whereby  the  specific  volume  of  the  solid  is  varied. 

5)  Determine  the  hot  spot  temperature  by  employing  Equations  (2.29)-(2.32). 
However,  the  hot  spot  temperature  is  only  calculated  after  complete  pore 
collapse  and  before  reaction  begins. 
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2.7  Numerical  Solution  Technique 

The  equations  presented  in  the  previous  sections  of  Chapter  2  completely 
define  the  hydrodynamics  and  thermodynamic  state  of  a  continuous,  porous 
medium.  A  finite  difference  technique  patterned  after  the  WONDY  code  [27]  was 
utilized  to  solve  the  system  of  equations.  At  the  initial  time,  t  =  0,  the 
bed  of  porous  propellant/explosive  is  discretized  into  J  cells  labeled  from 
left  to  right  as  j  =  1,2,3, ...,J.  The  thermodynamic  properties,  pressure, 
temperature,  internal  energy,  and  specific  volume,  are  assumed  to  be  constant 
over  the  width  of  each  cell.  At  the  boundaries  of  the  cells,  the  parameters, 
particle  velocity,  and  spatial  locations  are  defined.  Figure  2.4  illustrates 
the  location  of  the  points  of  evaluation  of  the  variables  in  time  and  spatial 
frames,  where  the  steps  in  time  are  represented  chronologically  as  n-1,  n, 
n+1.  Furthermore,  Table  [2.3]  presents  the  sequence  in  which  the  fluid  motion 
equations  are  solved,  which  are  written  in  their  finite  difference 
approximation  form.  Steps  (1-8),  in  Table  [2.3],  are  performed  to  model 
stress  wave  propagation  in  a  non  reactive  porous  material  while  steps  (1-4) 
and  (5+-9+)  outline  the  procedure  in  modeling  the  reaction  of  the  propellant. 

All  variables  are  initially  specified  for  each  cell,  located  throughout 
the  bed.  Following,  the  second  time  step  is  obtained  by  first  determining  the 
velocity  of  the  free  boundary  from  the  conservation  of  momentum  equation, 
written  as 


u 


n  +  V2 

V2  “ 


un  •  V2 

V2 


Pi)/”l 


(2.38) 


Time 


Eulerian  Coordinates 


Figure  2.4  Representation  of  finite-difference  cells  indicating  how 
the  variables  are  defined  with  respect  to  the  Lagrangian 
indices  and  Eulerian  coordinate  system.  Figure  taken 
from  Reference  [21]. 


Table  2.3 


Lagrangian  Finite-Difference  Equations 


u"  %-  Uj"  :l|*V4Ut"  ^  at"  -  -  (P*q)jn] 

(Mj  +  Vi> 


.  n  +  1 

hj  -  v. 


+  V2  _ 

tntl  -  t" 

V2-  hj 

_l^)/vj  =  constant 

n 

.  n  +  Vo  n  +  Vo 

j  -  v2+ 

st  2uj-i| 

n  +  1 

j  +  V2 

hj  -^/Hj 

where 

n+1 

qj(D 

=  0 

if  un 
J 

(ii) 

n+1 

qj(2)= 

ARV2 

n  +  1 
(hj  *  \  - 

where 

n+1 

=  0 

if 

n+1 

qj  (2 ) 

Vj  * 

( i  i  i ) 

n+l_ 

n+1 

n+1 

qj  ' 

qj(D 

+  qj (2 ) 

p"+1  =  f(a"+1)  determined  by  Equations  (2.16-2.20)  in  Section  (2.3) 

J  J 

.  n+1  n+1  ,  ,  »n+l 

where  ctj  3  /  (vs)j 
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6) 


Tn+^  determined  from  en+^  =  (e  )n+* 
J  J  s'j 


where  e"+1  =  e"  +  [(P>  q.j)n+1  +  (pj  +  9j)n]  ^7  +  aQ 


,  n+1 


.n+1 


i  n+1  • 


so 


<S>J  -  C*.(T  -  To>  *  ]  -  r  CvsToln  I.  ,  n+1- 

Vj 

0  n  n+1,  n  n+1, 

Av  .  ?,i  ,i  (,.i  -  ’i  1 

7  (vn  ♦  v"*1)2 


aQ  =  0  no  reaction 


7) 


,n+l. 


,  n+1 


Pn+1 


P'."1-  (P  )'.,T1=  Ire  (T'.,t1-  t  ) ( v  )n+1+  [ J ( ( v  )1T1)H/f(v  )nTY+  p 

j  s  j  1  vsv  j  o'v  s'j  av  L  1 ' j  7  J  r ' L  V v  s 7  j  J  0 


,n+l 


v  n+1 ,2, 


3) 


5)+ 


c.*  -  (v.ir1 1^ 


s  J 


3». 


[«Ps'i 


n+1 


partial  derivative  evaluated 
at  constant  temperature 


(i)  vT*  •  u”  -l/fcU'jUt"  *%*  At"  2  exp  [■ 


-T 


(T  )* 
H  j 


if  1  >  w"+1>  WQ 

—  J  —  b 


where  -  T  3  Equation  (2.35)  >_  -  ~ 

R 


(Tu)t  is  determined  from  Equation  (2.29) 

“  J 

(  ,  w"+1  =  w"-V2w"(Atn  +1/2+  Atn  ‘  1;2)  z  exp  [-  ^-f-\  if  WR  >  W*J+1  >  0 

J  J  J  R  T  B  J 

j 

where  CB1  T .  =  T? 

J  J 

see  explanation  on  page  45 
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®  V  ^  l<- r„ V  Tjnl  ,fTon<<TH>j 


Tn+^  determined  from  en+^=  Wn+^(e  )n+^  +  ( 1  -Wr?+ ^ )  ( e  )n+^ 
J  J  J  s'j  j  g  j 


w^re  (eg)f  -  cyg  (if1-  Tq) 


4q  .  (wj  -  Mj+1)  hdet 


Pn+1  =  R  Tn+1[l  +  0/(v  )n+1]/(v  )n+1 

j  j  L  g  j  J  g  j 


where  (vg)j+* 


n+1  n+1  ,  ,n+l 

li  ~M1  (,s>.i 

«-wJn+l) 


(P  )ntl 
s  J 


n+1 

Vj 

(»>"+1  (v0)"+1(«"+1-l)  V2 

If i-J +  3 - — J - 1  l 

^  n+lr*  n+1  *  JJ 

“j  Cs  «j  c9 


where  C  5  3  R  T1^  H  M  /  (v  )n+^ 

g  j  j  g  j 


*  The  localized  hot  spot  temperature  is  determined  after  complete 

compaction  and  the  values  used  for  specific  volume  correspond  to  the 


pressure  strength  of  the  shock  that  initiated  compaction. 
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with  q,  the  artificial  viscosity  which  will  be  described  later,  being  set 
equal  to  zero.  In  Equation  (2.38),  P^c  is  the  pressure  at  the  free  boundary 
condition  described  by  Equations  (2.1a  and  b).  Notice,  as  shown  in  Table 
[2.3],  the  conservation  of  mass  is  expressed  in  two  steps  (2  and  3),  thus 
enabling  a  solution  for  the  location  of  the  Lagrangian  spatial  coordinates, 
for  each  step  in  time,  explicitly.  Steps  (5-7)  and  (6+-8+)  have  to  be  solved 
implicitly,  therefore  an  iterative  technique  must  be  utilized  for  their 
solution,  while  the  remaining  steps  are  all  solved  explicitly. 

To  insure  stability  of  the  solution  over  time,  an  expression  for  the  time 
step  was  taken  from  Reference  [27]  and  is  written  here  as 


n+3/2_ 
4tj  • 


u  ,  n+1 
Ktl  Ah. 


VARB  +  { (VARB )  +  (Cn+1)  f2 

3 


VARB  =  B2C"+1  +  2  B \  Ahj+1 


(2.39a) 


(2.39b) 


where  K,  8-^  and  B2  are  constants  while 


„  ,  n  n+l> 

/.  I  JH  )  =  .  1 J.  1  -VJ _ \ _ n 

v  H  (vJ  +  vfht"  +V2 


(2.39c) 


Furthermore,  the  time  step  is  never  assigned  a  value  greater  than  five 
nanoseconds,  which  is  to  insure  stable  solution  during  decomposition. 

Also  to  insure  stability  of  solution,  q,  the  artificial  viscosity  term  was 
introduced  to  conservation  of  momentum  and  energy  equations  of  the  mechanical 
mixture. 

Artificial  viscosity  is  employed  in  the  numerical  finite-difference 

dP 

technique  to  prevent  an  infinite  slope  of  the  pressure  wave,  ,  which  would 
asymptotically  occur  when  a  shock  wave  is  formed.  Essentially,  by  utilizing 


v-v! 


.  \  1 
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artificial  viscosity,  the  shock  wave  is  spread  over  a  finite  number  of 
Lagrangian  cells.  A  time-intensive  effort  was  required  before  a  proper 
expression  for  artificial  viscosity  was  *ound  which  gave  reproducible 
results.  The  expression  for  artificial  viscosity,  presented  in  Table  [2.3]  as 
step  (4),  is  a  combination  of  a  quadratic  equation,  taken  from  Reference  [28] 
(which  essentially  spreads  the  shock  wave)  and  a  linear  equation,  taken  from 
WONDY  [27]  and  subsequently  modified  (which  dampens  the  oscillation  that  might 
occur  throughout  the  bed). 

Sound  speed  is  introduced  only  to  provide  mathematical  closure  for  the 
artificial  viscosity  expression  and  Equation  (2.39).  Steps  (8  and  9+)  in 
Table  [2.3]  represent  the  expressions  which  are  assumed  to  sufficiently 
describe  the  sound  speed  in  a  porous  and  two  phase  material,  respectively. 

The  numerical  solution  technique  described  above  was  incorporated  into  a 
code  given  in  Appendix  C. 
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CHAPTER  3 


NUMERICAL  RESULTS 

This  chapter  will  present  several  cases  which  capture  a  five-part 
scenerio  of  Deflagration  to  Shock  to  Detonation,  as  illustrated  in 
Figure  1.10.  All  the  cases  that  will  be  presented  numerically  model  a 
ten  centimeter  long  bed  of  HMX,  either  initially  homogeneous  or 
initially  porous.  The  bed  was  discretized  into  two  hundred  Lagrangian 
finite-difference  cells,  thus  corresponding  to  an  initial  grid  spacing 
of  five  tenths  of  a  millimeter,  (ax=0.05  cm).  Earlier  studies,  by 
Coyne,  Butler  and  Krier  [10],  employed  the  same  initial  grid  spacing  for 
modeling  stress  wave  formation  and  propagation  in  an  explosive  medium. 

In  addition,  the  values  of  the  coefficients  for  artificial  viscosity  are 
ARVl=2.0  and  ARV2=0.1,  which  are  used  for  all  cases  shown  in  the 
sections  to  follow  (unless  explicitly  stated  otherwise).  Also,  as 
stated  in  Chapter  2,  the  C82  combustion  model  is  utilized  unless 
explicitly  stated  otherwise.  All  propellant  properties  are  listed  in 
Appendix  A. 

The  following  sections  will  discuss  and  present  computed  results 
for  shock  initiation  to  detonation  for  both  homogenous  and  porous 
materials.  Validation  of  the  model,  based  on  several  typical  cases  of 
"ramp"  wave  initiation  to  detonation  in  a  porous  material,  is  also 
included.  In  addition,  a  comparison  of  the  two  combustion  models  (CB1 
and  CB2),  is  presented,  including  calculations  for  different  parameters. 


3.1  Shock  Initiation  of  Detonation 

The  first  calulations  considered  a  purely  homogeneous  bed  of  HMX.  Since 
the  bed  is  purely  homogeneous  (no  voids  exist),  the  characteristic  burn 
temperature  is  represented  simply  as  the  temperature  of  the  bed.  Also 
different  from  the  general  void  containing  cases,  analysis  for  a  homogeneous 
material  requires  an  artificial  viscosity  formulation  documented  in  Reference 
[271.  A  shock  strength  of  10  GPa,  i.e.,  P*=10  GPa  and  t*=0.02  ysec,  was  used 
to  initiate  the  stress  waves  that  produce  detonative  reactions  in  the  bed.  A 
different  criteria  was  needed  to  initiate  reaction  here,  since  the  temperature 
produced  by  shocking  the  homogeneous  medium  was  not  sufficient  in  initiating 
reaction  utilizing  the  explosion  model.  If  a  temperature  greater  than  450  K 
is  reached,  a  value  of  0.85  is  assigned  for  the  mass  fraction,  thus 
artificially  inducing  reaction.  Subsequently,  a  detonation  was  produced  in 
the  bed,  as  clearly  shown  in  Figure  3.1,  which  exhibits  the  pressure  profile 
of  the  detonation  wave  as  it  propagates  through  the  bed.  Listed  in  Table 
[3.11,  are  the  CJ  parameters  predicted  by  the  computer  code  (presented  in 
Appendix  C)  along  with  the  CJ  parameters  predicted  by  TIGER  [231.  Although 
the  CJ  parameters  predicted  by  code  do  not  agree  exactly  with  those  predicted 
by  TIGER,  the  results  are  acceptable. 

Problems  arose,  however,  when  trying  to  produce  a  shock  wave  in  a  porous 
material.  Figure  3.2,  taken  from  Reference  [21],  shows  the  velocity  magnitude 
of  the  free  boundary  as  a  function  of  time  for  a  homogeneous  case  (a0  =  1), 
and  two  porous  cases  (aQ  =1.2,  and  aQ  =1.4),  produced  by  a  "ramp"  wave  with 
P*  =  3  GPa  and  t*  =  40  usee.  One  can  clearly  see  a  "jump"  in  the  velocity 
magnitude  for  both  porous  cases  contrasting  with  the  homogeneous  case  which 
shows  a  smooth  transition  to  a  steady  value.  Coyne,  Butler  and  Krier  [10] 
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Figure  3.2  Velocity  of  the  left  boundary  as  a  function  of  time 

from  the  finite  difference  initial  porosities  compared 
with  the  velocity  of  a  one-dimensional  regressive  burn, 
r  =  0.001  mm/usec  P  (GPa).  Figure  taken  from  Ref.  [21]. 


attributed  the  jumps  in  velocity  magnitude  to  be  caused  by  the  pore  collapse 
model  and  not  by  numerical  inaccuracies.  Therefore  to  introduce  a  shock  wave 
into  a  porous  explosive,  a  specific  "ramp''  wave  is  required,  while  considering 
the  propellant  bed  to  be  inert  until  the  compression  wave  has  achieved  the 
desired  strength.  A  condition,  P*  =  ^  GPa  and  t*  =  10  psec,  was  used  to 
create  the  shock  wave  in  a  porous  explosive  case,  aQ  =  1.5322.  Figure  3.3 
illustrates  the  pressure  profiles  leading  to  a  steady  detonation.  The 
predicted  steady  CJ  values  for  pressure,  specific  volume  and  temperature  are 
19.49  GPa,  0.5544  cc/g  and  4597  K  with  a  corresponding  detonation  velocity  of 
7.082  mm/vjsec,  values  which  are  proper  for  HMX  with  a  density  of  pQ  =  1.24 
g/cc,  reduced  from  the  crystalline  density  of  p0  =  1.90  g/cc  when  a0  =  1. 

3.2  Validation  of  the  Model 

As  mentioned  in  Chapter  1,  accumulating  accurate  experimental  data  for 
shock  initiating  a  porous  explosive  is  difficult.  However,  J.  J.  Dick  [11] 
was  successful  in  obtaining  "Pop  plot"  data  through  extensive  runs  and 
statistical  analysis  for  samples  of  porous  HMX,  namely,  a0  =  1.5322.  Several 
calculations  were  made  here  assuming  a  porosity  identical  to  the  one  used  in 
Reference  [11]  by  J.  J.  Dick,  varying  the  shock  strength  in  order  to  produce  a 
theoretical  "Pop  plot".  For  comparison,  Dick's  results  are  displayed  in 
Figure  3.4  as  solid  lines,  along  with  the  computed  results  (shown  by  the 
dashed  line).  It  is  clear  from  Figure  3.4,  that  although  the  calculated 
run-up  to  detonation  distances  do  not  fall  in  the  uncertain  range  of  Dick's 
experimental  results,  the  model  does  predict  an  almost  identical  slope.  The 
quantitative  disagreement  could  stem  from  utilizing  calculated  decomposition 
times  for  the  hot  spots  based  on  experimental  data  obtained  for  a  similar  yet 
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4  A  comparison  of  numerically  predicted  and 
experimently  observed  [11]  run-up  distance 
versus  input  pressure  for  a  porous  bed  of 
HMX  (a  =1.5322) 


different  propellant,  HNS.  Futhermore,  Hayes'  hot  spot  model  was  based  on 
experimental  results  for  porosities  ranging  from  1.0875  <  aQ  <  1.2518,  where 
Dick's  experiment  and  the  calculations  were  conducted  for  a  higher  porosity, 
ct0  =  1.5322.  Finally,  a  true  shock  can  never  be  precisely  modeled  since  the 
stress  wave  must  be  spread  over  several  grid  elements,  a  discrepancy  in  the 
predicted  results.  Nevertheless,  within  the  limitations  of  the  model,  the 
results  predicted  are  exceptionally  good. 

3.3  Typical  Cases 

Several  cases  were  chosen  with  various  maximum  input  stresses, 
characteristic  rise  times  and  porosities  to  test  the  model  trends.  The  first 
case  considered  an  initial  porosity  of  o0  =  1.1176,  initiated  by  a  "ramp"  wave 
having  a  maximum  input  pressure  of  P*  =  2  GPa  and  characteristic  rise  time  of 
t*  =  10  psec.  Figures  3. 5-3. 9  illustrate  the  distribution  history  of  several 
key  parameters,  namely  P,  a,  (1-W),  T^,  and  T.  For  this  particular  case,  six 
specific  segments  in  time  were  chosen  for  illustration  purposes.  The  first 
time  shown  in  all  the  five  separate  figures  is  10  psec.  By  viewing  Figure 
3.5,  one  can  see  that  the  compression  wave  has  propagated  into  the  bed  to  a 
distance  of  1.5  cm.  Notice  even  though  the  left  boundary  has  reached  the 
maximum  input  stress,  shown  approximately  at  a  location  of  2  mm  into  the  bed, 
the  shock  front  has  not  fully  developed.  Figure  3.6  illustrates  the  closure 
of  the  voids  as  a  result  of  the  compression  wave,  or  in  other  words,  displays 
the  formation  of  the  solid  plug. 

One  can  see  from  Figure  3.7,  at  10  psec,  that  the  reaction  has  not 


commenced.  As  time  progresses  the  compression  waves  coalesce,  strengthening 
the  shock  front  which  in  turn  deposits  larger  amounts  of  (irreversible)  energy 
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into  the  voids  prompting  an  increase  in  the  hot  spot  temperature.  Figure  3.8 
exhibits  this  steady  increase  by  the  positive  slope  of  the  hot  spot 
temperature  profile.  Furthermore,  one  can  see  that  the  hot  spot  temperature 
is  much  greater  then  the  material  temperature  of  the  bed,  the  latter  depicted 
in  Figure  3.9. 

As  time  progresses  to  11  usee,  Figure  3.7  shows  reaction  has  begun.  This 
occurrence  can  also  be  viewed  by  an  increase  in  porosity  (alpha),  as  displayed 
in  Figure  3.6  at  that  instant.  Induced  by  the  the  initiation  of  propellant 
decomposition  an  increase  in  strength  of  the  compression  front  occurs,  as  seen 
in  Figure  3.5.  The  effect  of  the  strengthened  compression  front  is  a  rise  in 
the  hot  spot  temperatures  as  illustrated  in  Figure  3.8.  Also  because  of  the 
reaction,  there  is  a  predicted  increase  in  the  bulk  temperature,  exhibited  in 
Figure  3.9  at  t  =  11  ysec. 

At  the  third  time  presented,  1?  ysec,  the  degree  of  reaction  has 

increased,  again  shown  in  Figure  3.7  and  in  Figure  3.6.  With  the  increase  in 

the  degree  of  reaction,  the  chemical  energy  further  contributes  to  strengthen 
the  compression  front,  as  illustrated  in  Figure  3.5.  In  turn,  a  stronger 

shock  strength  prompts  a  higher  hot  spot  temperature,  depicted  in  Figure  3.8, 

and  the  additional  chemical  energy  also  stimulates  a  rise  in  the  material 
temperature,  as  shown  in  Figure  3.9. 

At  13  ysec,  the  decomposition  of  the  propellant  is  complete  in  a  small 
region  of  the  bed  near  x  =  2  cm,  as  pictured  in  Figure  3.7.  The  porosity 
distribution,  displayed  in  Figure  3.6  (at  13  ysec),  shows  compression  ahead  of 
the  decomposition  region.  The  effect  of  all  of  this  is  a  right  moving 
detonation  wave  and  a  left  or  rearward  moving  retonation  wave,  as  illustrated 
in  Figure  3.5.  One  may  notice  that  the  retonation  wave  has  a  higher  peak 


pressure  than  the  detonation  wave,  since  it  must  travel  through  a  highly 
compressed  material,  whereas  the  detonation  wave  must  propagate  through  a 
porous  material.  The  substantial  rise  in  the  pressure  front  subsequently 
causes  an  even  higher  hot  spot  temperature,  due  to  the  overabundant 
irreversible  energy  being  deposited  at  the  void  sites,  as  shown  in  Figure  3.8. 

At  the  fifth  time  shown,  t  =  14  psec,  the  region  of  complete 
decomposition  is  greatly  enlarged,  as  illustrated  in  both  Figure  3.6  and 
3.7.  The  chemical  energy  being  released  in  turn  prompts  the  propagation  of 
the  detonation  wave  further  into  the  propellant  bed,  as  seen  in  Figure  3.5. 
However,  the  retonative  reaction  has  ceased  propagating,  causing  a  steady 
decrease  in  the  left  moving  wave.  One  would  have  expected  that  the  retonation 
wave  would  have  propagated  all  the  way  to  the  left  boundary.  However,  the 
Arrhenius  kinetics  failed  to  sustain  the  retonation  wave,  just  as  the 
Arrhenius  kinetics  failed  to  initiate  reaction  in  homogeneous  material  shocked 
to  10  GPa.  Hayes  [17]  came  to  the  conclusion,  for  a  porous  material,  that 
scissing  of  molecular  chemical  bonds  induced  the  experimentally  observed 
decomposition  times.  It  is  therefore  postulated  that  a  shock  of  substantial 
strength  may  also  change  the  molecular  structure  of  a  voidless  propellant 
enhancing  reaction.  Also  at  time,  t  =  14  psec,  Figure  3.8  shows  that  the 
detonation  wave  produces  a  steady  hot  spot  temperature.  The  oscillations  seen 
in  the  hot  spot  temperature  are  a  result  of  numerical  integration 
inaccuracies. 

At  the  final  time  shown,  t  =  15  usee,  Figure  3.5  illustrates  the  steady 
detonation  wave  continuing  to  propagate  even  further  into  to  the  bed  at  a 
predicted  CJ  pressure  of  24.99  GPa  and  a  corresponding  CJ  temperature  of 
3923  K,  shown  in  Figure  3.9. 
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In  addition  to  the  figures  just  presented,  a  locus  of  stress  and  reaction 
fronts  are  shown  in  Figure  3.10.  The  dashed  line  depicts  stress  wave 
propagation  into  the  bed,  initiating  the  detonation  t  =  12.41  psec,  at 
x  =  2.05  cm.  The  solid  lines  represent  the  locus  of  the  right  and  left 
traveling  detonation  fronts.  Figure  3.10  also  shows  the  termination  point  of 
the  retonation  wave.  An  apparent  change  in  velocity  may  be  seen  in  the  left 
traveling  wave  at  the  termination  point.  The  slope  of  the  solid  lines 
correspond  to  a  detonation  and  retonation  velocity  of  D  =  8.757  mm/psec  and 
R  =  9.135  mm/psec. 

A  second  case  is  now  presented  for  a  material  with  a  higher  initial 
porosity,  namely  aQ  =  1.1875,  and  with  an  increased  strength  "ramp"  wave 
P*  =  4  GPa  and  t*  =  10  psec.  Figures  3.11-3.15  depict  the  pressure,  porosity, 
mass  fraction  reacted,  hot  spot  temperature,  and  bulk  temperature  distribution 
for  fixed  time  segments.  One  can  see  from  Figure  3.11  that  the  detonation 
occurs  just  behind  the  compression  front  and  in  the  following  time  the  stress 
front  is  overtaken  by  the  detonation  wave. 

Different  from  the  previous  case  presented,  no  retonation  takes  place. 
However,  a  left  moving  compression  can  be  seen  in  Figure  3.11.  If  a  chemical 
change  does  occur  when  a  homogeneous  material  is  strongly  shocked  (enhancing 
the  sensitivity  to  reaction),  as  postulated  earlier,  a  retonation  would  have 
been  caused.  Also  illustrated  in  Figure  3.11  is  the  steady  state  detonation 
wave  having  a  CJ  pressure  of  24.39  GPa  and  a  corresponding  CJ  temperature  of 
4173  K  depicted  in  Figure  3.15.  The  physical  plane,  shown  in  Figure  3.16, 
illustrates  a  detonation  velocity,  D  =  8.255  mm/psec,  but  no  retonation 
wave.  The  detonation  occurs  at  a  run-up  distance  of  x  =  4.83  mm  at 
t  =  4.96  psec. 
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Predicted  unreacted  mass  fraction  distribution  history  for  a  porous  bed  of  HMX 
(a  =1.1875,  P*=4  GPa,  t*=10  psec). 
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The  third  case  treats  a  material  with  a  relatively  high  initial  porosity. 


namely  oQ  =  1.4615  but  a  relatively  weak  ramp  input  wave,  P*  =  2  GPa  and 

t*  =  20  ysec.  Again  the  five  parameters  P,  a,  ( 1 -W ) ,  T^,  and  T  are 

illustrated  in  Figures  3.17-3.21,  respectively,  as  time  progresses.  Similar 
to  the  previous  case  (a0  =  1.1875),  no  retonation  is  predicted  to  occur,  but  a 
left  traveling  compression  can  be  seen.  In  this  particular  case  a  steady  CJ 
pressure  of  21.01  GPa  and  CJ  temperature  of  4616  K  are  predicted.  One  can  see 
a  peak  in  the  hot  spot  temperature  in  the  regime  of  initial  detonation.  Since 
a  constant  hot  spot  temperature  would  be  expected  after  a  detonation  wave  is 

initiated,  the  predicted  peak  in  the  hot  spot  temperature  (Figure  3.20)  is 

assumed  to  be  caused  by  inaccuracies  in  the  numerical  steps  used  to  obtain  a 
hot  spot  temperature.  The  physical  plane,  shown  in  Figure  3.22,  again 
presents  the  location  and  time  of  the  occurrence  of  detonation,  x  =  2.68  cm 
and  t  =  21.55  ysec.  The  slope  of  the  solid  line  represents  the  detonation 
velocity,  0  =  7.24  rnm/ysec. 

3.4  Numerical  Accuracy  Test 

In  addition  to  the  cases  just  shown,  a  comparison  was  made  of  the 
pressure  profile  for  two  different  initial  discretized  cell  dimensions, 
depicted  in  Figure  3.23,  for  HMX  with  an  initial  porosity,  a0  =1.2667,  and  the 
imposed  stress  input  condition,  P*  =  5  GPa  and  t*  =  10  ysec.  One  can  see  that 
the  location  of  the  compression  wave  for  both  initial  grid  spacings  is 
approximately  the  same  for  the  first  two  steps  illustrated.  However,  the 
difference  in  the  pressure  magnitudes  for  the  two  initial  grid  spacing  is 
sufficient  to  cause  a  detonation  wave  to  be  predicted  later  in  the  larger 
initial  grid  spacing.  The  corresponding  location  and  time  of  detonation  for 
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the  initial  cell  dimensions,  ax  =  0.1  cm  and  ax  =  0.025  cm,  is  x  =  5.189  mm, 
t  =  10.65  usee,  and  x  =  5.668  mm,  and  t  =  10.59  psec,  respectively.  Also 
notice  at  t  =  14  psec,  the  larger  initial  grid  spacing  does  not  capture  the 
peak  pressure. 

For  the  majority  of  cases  studied  here  (and  due  to  a  limited  amount  of 
computer  funds)  a  minimum  initial  cell  size  was  chosen  as  Ax  =  0.05  cm.  This 
value  allowed  us  to  retain  essential  features  of  the  shock  initiation 
physics.  Furthermore,  for  the  parametric  studies  carried  out  here  the  basic 
trends  will  still  be  accurate.  Nevertheless,  quantitative  results,  such  as 
the  point  of  detonation  initiation  may  be  somewhat  in  error. 


3.5  Comparison  of  the  Combustion  Models 

Shown  in  Figure  3.24  is  a  comparison  of  the  predicted  shock  run-up  to 
detonation  length  for  combustion  models  CB1  and  CB2,  for  an  input  condition  of 
P*  =  2  GPa  and  t*  =  10  psec.  The  two  models  agree  exactly  for  the  last  four 
high  alpha  values  studied.  However,  there  is  a  considerable  difference  for 
the  low  porosity  case,  a0  =  1.1176.  In  addition,  utilizing  combustion  model 
CB1,  a  detonation  was  not  produced  for  the  low  initial  porosity  of 
a0  =  1.0566,  when  P*  =  2  GPa  and  t*  =  10  nsec. 

Figure  3.25  illustrates  the  predicted  pressure  profile  for  an 
a0  =  1.1176,  employing  CB1.  Notice  that  the  detonation  is  predicted  to  occur 
at  a  much  later  time,  and  a  substantial  distance  behind  the  compression  front, 
as  compared  to  the  same  case  employing  CB2,  previously  shown  in  Figure  3.5. 
However,  the  essential  features,  that  is,  detonation  and  retonation  velocities 
with  the  same  CJ  properties  are  predicted.  Figure  3.26  shows  the  predicted 
temperature  profile  with  the  resulting  CJ  conditions,  while  the  physical 
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Predicted  temperature  distribution  history  for  a  porous  bed  of  HMX 
utilizing  combustion  model  CB1  (a  =1.1176,  P*=2  GPa,  t*=10  psec). 


plane,  presented  in  Figure  3.27,  shows  the  detonation  and  retonation 
velocities. 

Since  CB2  was  able  to  produce  detonations  in  both  low  porosity  cases 
studied,  and  only  slight  disagreements  were  found  in  the  predicted  CJ 
properties  between  the  models,  combustion  model  CB2  is  the  logical  choice. 
Table  3.2  lists  the  CJ  parameters  along  with  the  detonation  velocities  for  all 
the  porosities  studied.  Comparing  Tables  [2.1]  with  Table  [3.2],  one  can  see 
that,  although  agreement  between  the  TIGER  equilibrium  calculations  [23]  and 
the  predicted  results  presented  in  Table  [3.2]  are  not  exact,  the  same 
qualitative  trends  are  present.  However,  problems  did  occur  in  predicting  the 
proper  Pqj.  The  proper  P^j  can  be  predicted,  with  some  effort,  by  readjusting 
the  artificial  viscosity  coefficients  for  each  porosity. 

Furthermore,  from  the  results  presented  in  Figure  3.24,  one  can  clearly 
conclude  that  alpha  is  not  a  dominant  factor  in  the  effective  run-up  distance 
to  detonation  for  the  interval  1.15  <  a0  <  1.4615.  One  would  expect  a 
decrease  in  run-up  distance  to  detonation  as  the  porosity  increases.  However, 
for  a0  =  1.4615  an  actual  increase  in  run-up  distance  is  predicted,  although 
the  time  to  detonation,  displayed  in  Figure  3.28,  always  decreases 
monotonical ly  for  increasing  porosity.  A  possible  cause  for  the  increase  in 
run-up  distance  can  be  attributed  to  the  amount  of  initial  void  volume. 

Figure  3.29  portrays  the  location  of  the  left  boundary  as  a  function  of  time, 
utilizing  an  input  condition  with  P*  =  3  GPa  and  t*  =  40  ysec.  One  can 
clearly  see,  an  effect  of  the  initial  void  volume,  namely  that  for  the  greater 
porosities  the  left  boundary  shifts  further  into  the  bed.  Since  the  run-up 
distance  is  measured  from  the  initial  location  of  the  left  boundary  to  the 
location  of  the  first  Lagrangian  finite  difference  cell  to  detonate,  the 


0*0T 


retonation  fronts  (a  =1.1176,  P*=2  GPa,  t*-10  psec). 
his  particular  case  Bombustion  model  CB1  was  utilized. 
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TABLE  3.2 

■  PREDICTED  CJ  PARAMETERS 


E 

B 


a 


1.0566 

1.1176 

1.1875 

1.2667 

1.3571 

1.4615 


pTo 

(g/cc) 


1.8 

1.7 

1.6 

1.5 

1.4 

1.3 


rCJ 

(GPa) 


23.76 

24.99 

24.03 

24.13 

21.37 

21.16 


'CJ 

(K) 


3532 

3923 

4151 

4371 

4430 

4629 


VCJ 

(cc/g) 


0.4571 

0.4752 

0.4873 

0.5006 

0.5206 

0.5324 


D 

(mm/ysec) 

10.48 
8.759 
8-.  480 
8.133 
7.685 
7.177 


*  All  results  produced  with  a  P*  =  2  GPa  and  t*  -  10  ysec. 


Time  to  Detonation  (cm) 


10.0  L_ 

1.0 


P*=2.0  GPa 
t*=10  Msec 


■  CB2 
- CB2 


1.2  1.3 

a  (POROSITY) 


Figure  3.28  Predicted  time  to  detonation  versus  porosity 
utilizing  combustion  model  CE2  (P*=2.0  GPa, 
t*=10  usee) . 
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100. 


Position  of  Left  Boundary,  x0(t) 

Pressure  Boundary  Condition 
0<[P=  (0.075  GPa/fjLs)  t]<  3  GPa 


:ii] 


TIME  C  MICRQSEC .  ) 

Figure  3.29  Position  of  the  left  boundary  as  a  function 

of  time  from  the  finite  difference  calculations 
for  three  different  initial  porosities,  taken 
from  Reference  [21].  (P*=3  GPa,  t*=40  ysec). 


compaction  occurring  in  the  compression  front  ahead  of  the  point  of  detonation 
could  be  sufficiently  large  resulting  in  an  increased  run-up  distance 
calculated  for  large  porosity  conditions.  Another  probable  cause  may  have 
been  the  over  extension  of  Hayes'  hot  spot  theory  [17],  since  the  highest 
porosity  he  studied  was  an  aQ  «  1.252. 

3.6  Parametric  Studies 

To  conclude  the  research,  two  parametric  studies  were  made.  The  first 
compared  the  run-up  distances  for  a  fixed  porosity  (a  -  1.2267)  to  various 
stress  input  conditions,  while  the  second  compared  run-up  distances  for 
various  input  pressures  with  a  fixed  characteristic  rise  time  (t*  =  10  ysec) 
for  different  porosities.  Figure  3.30  presents  the  first  parametric  study. 

One  would  expect  a  greater  run-up  distance  to  correspond  to  a  longer 
characteristic  rise  time  for  a  specific  input  pressure,  since  the  compression 
front  will  reach  the  critical  pressure  needed  to  initiate  detonation  later  for 
the  longer  characteristic  rise  time.  Figure  3.30  demonstrates  the  occurrence 
of  a  longer  run-up  distance  for  a  specific  input  pressure  corresponding  to  a 
longer  characteristic  rise  time.  The  "Pop  plot"  (Figure  3.4)  shows  that  for 
shock  initiation  of  detonation  the  run-up  distance  increases  with  weaker  input 
pressures.  Similarly,  depicted  in  Figure  3.30,  for  "ramp"  wave  initiation  of 
detonation  the  run-up  distance  increases  with  decreasing  input  pressures. 

By  extrapolating  data  from  Figure  3.30  with  input  conditions  having  a 
P*  *  2  GPa,  a  long  run-up  distance  to  detonation  is  expected  for  slow 
characteristic  rise  times,  displayed  in  Figure  3.31.  Therefore,  for  low  input 
pressures  and  critically  long  characteristic  rise  times,  detonation  is  not 
expected  to  occur  in  a  ten  centimeter  bed,  eliminating  the  hazard  of  DSDT. 
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Predicted  run-up  distance- versus  input  press 
several  characteristic  rise  times  (a  =1. 2667) 
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Figure  3.31  An  extrapolation  to  determine  the  expected 
run-up  distance  for  relatively  long 
characteristic  rise  times  (a  =1.2667, 

P*=2  GPa) .  0 


Moreover,  Figure  3.32  shows  an  increase  in  run-up  distance  to  detonation 
to  correspondingly  weaker  input  pressures  for  four  specific  porosities, 
a0  ■  1.1875,  a0  =  1.2667,  a0  =  1.3571,  and  <*0  =  1.4615.  Evidence  of  numerical 
integration  errors  were  encountered  for  low  porosities,  aQ  <  1.15.  Therefore 
no  results  are  shown.  One  would  expect  a  higher  porosity  bed  to  be  more 
sensitive  to  detonation  than  a  lower  porosity  bed,  thus  resulting  in  shorter 
run-up  distances.  Therefore  the  curves  for  different  alphas  should  not  cross 
each  other.  However,  the  crossing  may  be  attributed  to  a  course  initial  grid 
spacing.  Although  there  exists  some  quantitative  errors  in  Figure  3.32,  one 
may  conclude  that  for  porosities,  in  the  interval  of  1.15  <  a0  <  1.5,  do  not 
take  a  dominant  role  in  the  run-up  distances.  Furthermore,  at  relatively  weak 
inputs,  P*  <  2  GPa,  the  run-up  distances  asymptote  to  high  values,  indicating 
that  no  detonation  would  be  predicted  for  input  pressures  lower  than  1.5  GPa. 


several  porosities  (t*=10  Msec). 


CHAPTER  4 


CONCLUDING  REMARKS 

The  previous  chapter  presented  several  numerically  produced  results  which 
depicted  the  five-part  scenario  of  Deflagration  to  Shock  to  Detonation 
Transition.  Although  some  integration  inaccuracies  in  the  predicted  results 
were  encountered,  the  results  clearly  showed  that  "Ramp"  wave  to  Detonation 
Transition  is  a  prominent  hazard  associated  with  porous  explosives.  One  of 
the  purposes  of  this  last  chapter  is  to  recommend  necessary  improvements  for 
further  utilization  of  the  code. 

4.1  Necessary  Improvements 

Two  areas  clearly  in  need  of  improvements  are:  (1)  the  need  for  a  better 
integration  scheme  to  define  the  shock  waves,  and  (2)  a  better  data  base  to 
calculate  the  hot  spot  temperature  and  subsequent  reaction  rates.  One  of  the 
first  steps  that  should  be  taken  to  accurately  define  the  shock  would  be  to 
reduce  the  amount  of  artificial  viscosity.  By  doing  this,  the  shock  wave  will 
be  spread  over  a  smaller  amount  of  grid  spaces.  If  additional  increased 
computer  funding  becomes  available,  it  is  recommended  that  the  initial  finite 
difference  cell  sizes  be  reduced  by  V? to  V4  .  The  result  of  both  reducing  the 
artificial  viscosity  and  initial  cell  size  will  be  better  representation  of 
the  shock  structure,  which  should  alleviate  most  of  the  numerically 
encountered  problems.  Furthermore,  a  more  careful  definition  of  the  shock 
structure  will  allow  a  more  accurate  calculation  of  the  hot  spot  temperature, 
since  one  of  the  errors  associated  with  the  determination  of  a  hot  spot 
temperature  stems  from  the  relatively  poor  numerical  representation  of  the 
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shock  structure.  Basically,  an  inaccurately  distributed  shock  wave  results  in 
an  error  in  the  evaluation  of  the  irreversible  energy,  since  the  compression 
process  is  performed  over  several  cells. 

Moreover,  since  the  reaction  of  a  porous  material  initiates  and  is 
controlled  in  the  initial  stages  of  decomposition  by  the  hot  spot,  a  better 
data  base  is  needed  to  accurately  determine  the  hot  spot  temperature  and 
decomposition  rates.  Although  the  Hayes  hot  spot  theory  [17],  incorporated  in 
the  code,  modeled  initiation  and  decomposition  exceptionally  well  for 
materials  with  high  porosities,  the  theory  began  to  collapse  when  alpha  was 
less  than  1.15.  Therefore,  it  is  suggested  that  research  should  be  conducted 
on  tne  formulation  of  a  model  for  the  localized  hot  spots  on  a  microscale, 
both  experimentally  and  theoretically.  Also,  directly  related  to  the' hot  spot 
temperature,  experimentally  observed  decomposition  times  should  be  measured 
for  a  wider  spectrum  of  porosities.  With  a  better  understanding  of  the 
formation  and  decomposition  of  hot  spots,  an  improvement  in  the  quantitative 
results  should  be  possible. 

In  addition  to  the  two  areas  just  mentioned  in  need  of  improvement,  a 
dynamic  pore  collapse  model  outlined  in  Reference  [20]  should  be  implemented 
in  the  code  either  to  validate  or  contradict  the  use  of  the  static  pore 
collapse  model  now  employed.  Although  Kooker  and  Anderson  [9]  found  in  their 
studies  that  the  static  pore  collapse  theory  sufficiently  modeled  compaction, 
the  inertial  and  viscous  terms  contained  in  the  dynamic  pore  collapse  model 
may  be  of  significance  in  defining  the  rapid  compression  arising  from  a  steady 
state  detonation  wave. 
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Although  the  code  did  not  successfully  produce  accurate  quantitative 
results  for  the  lowest  porosities  studied  here,  the  key  elements  of  the  five- 
part  scenario  were  modeled  predicting  a  steady  state  detonation  and,  in  some 
instances,  the  occurrence  of  a  retonation.  The  model  was  validated  by  a 
comparison  with  an  experimental ly  produced  run-up  length  versus  peak  shock 
variation,  the  so-called  "Pop  plot".  This  comparative  study  showed  the 
analysis  to  be  somewhat  conservative,  predicting  shorter  run-up  distances  for 
a  given  peak  pressure. 

However,  the  slope,  £SDy  vs  P*  was  matched.  A  parametric  study,  which 
varied  the  input  conditions  for  a  specific  porosity,  showed  that  longer 
characteristic  rise-time  resulted  in  a  longer  detonation  run-up  distance.  A 
second  parametric  study,  which  compared  run-up  distances  for  different 
porosity  explosives,  for  various  pressure  inputs  with  the  same  characteristic 
rise  time,  illustrated  that  porosity,  in  the  interval  1.15  <  a0  <  1.5,  was  not 
a  dominant  factor  in  the  effective  run-up  distance. 

The  granulated  bed/cast  explosive  configuration  is  representative  of  a 
rocket  motor  which  has  partially  fragmented.  Even  though  the  length  of  the 
fragmented  propellant  is  not  long  enough  to  detonate  from  an  accelerated 
convective  burn,  the  intact  propellant  may  shock  initiate  from  the  rapid 
pressure  rise  rate.  It  is  evident  that  "ramp"  wave  initiation  of  detonation 
is  truly  a  serious  hazard  to  contend  with,  demonstrated  by  the  short  run-up 
distances  corresponding  to  rapid  rise  rates  for  several  porous  cases. 

However,  by  the  extrapolation  of  data  from  the  first  parametric  study,  a 
substantially  long  run-up  length,  outside  the  dimensions  of  the  propellant 
bed,  would  result  for  relatively  long  rise-time  conditions. 


If  the  distance  needed  to  shock  initiate  the  cast  material  is  greater 
than  any  dimension  of  the  rocket  motor,  OSDT  is  obviously  impossible.  If  the 
solid  rocket  propellant  is  unavoidably  frangible,  then  future  formulations 
should  be  devised  so  that,  when  damaged,  the  propellant  breaks  into  larger 
fragments  (on  the  order  of  V2  millimeter).  Flame  spreading  and  convective 
burning  in  such  effectively  large  fragments  (but  smaller  surface-to-volume) 
would  result  in  a  local  dP/dt  that  would  produce  larger  t*. 

As  a  final  conclusion,  one  can  clearly  state  that  an  alternative 
methodology  for  transition  to  detonation  other  than  the  direct  acceleration  of 
the  convective  burn  front,  is  a  Ramp  Initiated  Shock  to  Detonation  Transition 
in  an  impermeable  but  porous  propellant. 


APPENDIX  A 


HMX  PROPERTIES  [9,22,23,26,29] 


Explosive:  0CTAHYDR0-1,3,5,7-TETRANITR0-1,3,5,7-TETRAZ0CINE 


Formula:  C^HgNgOg 
Initial  Homogeneous  Specific  Volume 
Ambient  Homogeneous  Sound  Velocity 
Initial  Yield  Stress 
Initial  Shear  Modulus 
Gruneisen  Coefficient 
Activation  Temperature 
Frequency  Factor 
Product  Gas  Constant 
Specific  Heat  (constant  volume) 
of  the  solid 

*Specific  Heat  (constant  volume) 
of  the  product  gas 

*Heat  of  Oetonation  HDET=  ^*91 

*Detonation  Velocity 


v$0  =  0.5263  cc/g 
cSQ  =  0.2642  cm/ysec 
Y0  =  51.7  MPa 
Gq  =  3.516  GPa 
r  =  1.1 

—  =  14400  K 

TT 

Z  =  6.9  x  1010 

R  =  2870740  £££ 
gK 

Cvs  =  1.5  J/gK 

CvG  =  [2.4-0.28(  -  1.3)]-~ 

To  y 

4.33(-i -  1.3)2-  0.934(— ^ - 1.3)]^ 

vTo  vTo  g 

0  =  3 . 64 (— -^ - 1.3)  +  6.98  mm/ysec 

vTo 


Covolune  Correction  Term 


6  is  listed  in  Table  [2.3] 
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Nonlinear  volume-dependent  function  obtained 
from  shock  Hugoniot  experiments  [22]  c 


J(x)  =  7.57x2  +  13.33x3  +  18.04x4 
+  2.828x5  +  24.01x6  +  278. 3x7 
+  383. 6x8 


where  x  =  — - - 1 

s 


*Data  fits  of  CJ  data  predicted  by  TIGER  [23]. 
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APPENDIX  B 


RECIPROCITY  RELATIONS 

Since  very  few  mechanical  engineers  utilize  the  concept  of  Helmholtz  free 
energy,  a  step  by  step  derivation  for  pressure  and  specific  internal  energy 
from  the  free  energy  definition  will  be  presented.  Helmholtz  free  energy  is 
defined  to  be 

4>  =  e  -  Ts  (8.1) 

where  4*  represents  Helmholtz  free  energy,  e,  specific  internal  energy,  T, 
temperature,  and  s,  specific  entropy.  By  taking  the  derivative  of  Equation 
(B.l)  the  following  expression  is  obtained 

d4»  =  de  -  Tds  -  sdT  (8.2) 

Making  use  of  an  important  thermodynamic  relation 

de  =  Tds  -  Pdv  (8.3) 


where  P  represents  pressure  and  v  specific  volume,  a  substitution  can  be  made 
in  Equation  (8.1)  yielding 


Evaluating  Equation  (B.4)  at  constant  temperature  an  expression  for  pressure 
in  terms  of  free  energy  is  arrived  at,  i.e. 


In  a  similar  manner,  an  expression  for  specific  internal  energy  can  be  derived 
by  first  evaluating  Equation  (B.4)  at  constant  specific  volume,  bringing  forth 

cffJv  ■  -s  <B-6> 

Following  a  rearrangement  of  the  definition  of  free  energy  and  then  using 
Equation  (B.6),  specific  internal  energy  can  be  expressed  as 

«  ■  *  -  T  (f*)„  (B.7) 

Although  the  derivations  may  seem  somewhat  trivial,  including  the  fundamentals 
here  may  assist  in  the  interpretation  of  the  model  and  the  results  for  any 
given  explosive.  See  Reference  [30]  for  additional  information. 
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